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Abstract. We study the Potts model on locally tree- like random graphs of arbitrary 
degree distribution. Using a population dynamics algorithm we numerically solve 
the problem exactly. We confirm our results with simulations. Comparisons with 
a previous approach are made, showing where its assumption of uniform local fields 
breaks down for networks with nodes of low degree. 
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Recent years have seen intensive interest in the statistical physics community in 
complex networks. Random graphs with a given degree distribution (often also called 
the configuration model) [H E], small world graphs |S1 IH E] and various scale-free 
graphs based on preferential attachment [HI EI have been extensively studied |H1 1^ ITU]. 
Many researchers have also considered processes taking place on complex networks, 
such as percolation [TT[ [T21 UHl , disease spreading [31[Tlj, self-organised critical models 
[T31 Uni HZ], and various games [IHl d]. These processes have in general been studied 
in the past on regular lattices, the current interest being in how the topology of the 
complex network affects the processes. The Ising model and its generalisation the q- 
state Potts model have also recently been studied on random graphs of a given degree 
distribution [201 HH 122] • These models are standard models of statistical physics, 
displaying continuous (for Ising) and first order (g > 2) phase transitions [22] in mean 
field. Their behaviour on complex networks, as well as being of interest in themselves, 
also gives clues as to how more complicated processes will behave. 

The solution of the Potts model on a Bethe Lattice has been known for some time 
[21]. More recently, Dorogovtsev, Goltsev, and Mendes have extended the solution to 
random trees [21] . This includes the configuration model for large network sizes since in 
these cases the graph is locally tree-like. Their approach, which is tailored to scale-free 
random graphs, relies on an effective field approximation which accounts for the effect 
of the nodes with the highest degree. This is indeed able to capture the correct nature 
of the phase transition and provides a qualitatively correct behaviour. 

In this paper we describe a method of numerically solving for the local fields 
without making any approximation. This shows that, at odds with the effective field 
approximation [21], local fields and the magnetization indeed depend on the local 
geometry and hence fiuctuate with it. The effects show their relevance, for example, in a 
recent study where the Potts model on a random graph was used to model coordination 
games on an evolving network |25]. 

We compare our results with the method of [21] and with simulations. We find 
that, as expected, the approximation of [21] gives poor results for graphs containing 
nodes of low connectivity but works well when all nodes have large connectivity. 

1. Potts Solution 

The g-state Potts model has the Hamiltonian 



where the as take one of the integer values from 1 to g and the sum is over 

neighbouring nodes. Here the applied magnetic field b acts on the (arbitrarily chosen) 
state 1. The partition function is then: 



(1) 



<i,j> i 




(2) 
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where K = pj and H = ph, (3 = l/{kBT), kB, T being Boltzmann's constant and 
temperature respectively. The local magnetisation rrii is given by 

q- 1 

where the angled brackets denote the Gibbs average. The total magnetisation is then 
the average of over all nodes. 

We are interested in the Potts model on a (connected) random graph of n nodes 
with specified degree distribution P{k). These are locally tree-like graphs in the sense 
that the sub-network of all nodes at a distance less than i steps from any particular 
root node is very similar to a tree, for small i. In particular, this is true for i <^ d, 
where d is the diameter of the network, which is known 26 to behave as ~ logn if 
the second moment of P{k) is finite and d ~ log logn if the second moment diverges, 
but the first is finite. The strategy of the solution relies on the following idea: We take 
a particular site as the root and show that its statistical properties can be derived from 
those of their neighbors. In turn, the properties of the neighbors are defined by those of 
second neighbors and so on. As long as the graph is tree-like, i.e. for distances r -C o? 
from the root, the statistical properties of the sites involved in the r**^ iteration - those 
at distance r from the root - are independent, and are captured by a simple transfer 
matrix. If this transfer matrix has a fixed point, the statistical properties of bulk nodes 
will be described by the fixed point irrespective of the boundary conditions which are 
set at the leaves of the tree (i.e. at distance r ~ logn). This will be true, in particular, if 
the boundary conditions are drawn from the fixed point distribution which describes the 
statistics of bulk sites. Notice that this argument requires the diameter over which the 
network is tree-like be much larger than the correlation length over which the transfer 
matrix converges. Hence for finite graphs we expect the local tree approximation to 
perform poorly for graphs with d ~ log logn and also close to the critical point where 
the convergence of the iterative process is slowed down. 

For a tree, following the notations of Ref. [21], we can solve Z iteratively, we define 

Where the sum over / and m is only over nodes in the sub-tree whose root node is i, 
including i. Thus 

^ = ;^e^^-^n^?i,(ao) (5) 

o-Q i=l 

and 

/ V ^ exp(g) - Hill Xl,^ , . 

^"^'^ exp{H) + {q-l)UtiX,, 

where Xi^i = gi^i{a) / gi^i{l) . Likewise, we can define the quantity (7r,j(o"r-i) for any node 
i at a distance r from the root, as the partial sum of the partition function on the 
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sub-tree ensuing from that node, where (T^-i is the value of the spin at distance r — 1 
where the sub-tree is connected. This satisfies the recursion relation 

9rA^r-l) = J^e^^-— n 9r+lA^r) (7) 

CTr j = l 

where kr,i is the degree of node (r, i). Notice that Qr^i depends on the value of the spin 
(Jr-i at the previous level. Writing down these equations for ar-i = 1 and aj—i = a > 1 
and taking the ratio one finds an iterative equation for 



which reads 



a > 1 



(8) 



(9) 



In principle Xr,i should depend on a in Eq. (jH)). However it is easy to see that Xr,i is 
independent of a if the node (r, i) is a leaf of the tree. Because all Xr^i can be computed 
recursively with Eq. (p)) from the leaves back to the root then Xr,i is independent of 
a > 1 for all nodes. 

From here, making the change of variables hrj = — In(xrj), Dorogovtsev et. al. 
fn\ noted that the only dependence on h in the right hand sides of equations was 
from a sum over the independent hiS. They thus made the approximation 

k 

Y,K+i,,^k{K+i) + 0{Vk). (10) 

i=i 

It is possible instead to solve for the distribution of h values in the steady state. 
Let pr{h\k) be the distribution density of hr^i with kj-^i = k. This distribution satisfies 



Pr{h\k) 



where 



and 



Y(s) 



In 



k-l 

Yl dhjpr+i{hj)6 

^H+K + (g _ i)e-^ 
+ (g - 2 + e^)e- 



fe-i 

h-Y\Y.h, 



Pr{h)=Y.P{k)pr{h\k) 



:iii 



(12) 



(13) 



k=l 



is the distribution of h on the neighbours of a node, P{k) = kP{k)/ (k) being the 
distribution of degrees of neighbours of a node. Averaging over P{k) we find a recursion 
relation for the distribution of fields at the neighbour of a node, Pr{h) 



Pr{h) = Y.P{k) / l[dh,Pr+l{h,)S 

k=l ■^^°°j=l 



h-Y{Y.h, 



(14) 



Well inside the bulk of the tree, we expect that pr ^ p independent of r. Once p{h) is 
found, the distribution of local fields on any site with degree k is given by Eq. (fTTj). 
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2. Algorithm for Solution 

A simple way to solve equation (fT^ numerically is the following: Start from a population 
hi, i = l,...,MofM ^ 1 values of h. Evolve the population by iteration with the 
following procedure: 

• Draw at random k from the distribution P{k) 

• Draw k — 1 values of h at random from the population {hi} and sum them to get 

hsum- 

• replace a random member of the population by 

Iterate until convergence. The magnetisation is then found by averaging the local 
magnetisations on many nodes, 

• draw k from P{k) 

• draw k values of h from p{h) and sum them to get hgum and insert this into 

(mo) - exp(ij') - e^Y>{-hsum) ^^^-^ 

Iterate many times and average the results to get M{T). 

By contrast, the approximation of Dorogovtsev et. al. (21] gives: 



{h) =j:pik)Yi{k-i){h)) (17) 



k=l 

oo 



To find the general M{T) these two equations are also solved numerically, albeit for 
only one value of h rather than a whole population. 



3. Results and comparison with Simulations 

In the following we compare the results of our method with that of Dorogovtsev et. al. 
and with simulations. In particular, for zero applied field {H = 0), we study the g = 10 
state Potts model on random graphs with degree distributions P{k) that are bimodal 
or truncated power law, 

p{k) = ^-hl±^ (19) 

P{k) = Ak for kmin ^ k ^ kmax (20) 

where A is the normalisation. 

The random graphs for the simulations are produced in the following way: n nodes 
i are given degrees ki drawn from P{k), if J2 ^i is odd then a randomly chosen node has 
its degree increased by 1. Nodes of degree k are imagined to have k unconnected ends. 
Next, pairs of ends are chosen at random and linked iff they belong to different nodes 
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Figure 1. Plot of M{T) for the g = 10 Potts model on a random graph with bhiiodal 
degree distribution (see equation (|19|l l with J/ks — 1- Crosses are simulations for 
2 realisations of a network with n — 20, 000 nodes, starting from an initially fully 
magnetised state. Lines are the results of the 2 methods described, the method 
introduced here agreeing with simulations to high accuracy. The vertical line at T « 1.1 
is the transition from the metastable unmagnetised state to the magnetised state, both 
methods giving the same value of Tmetastabie ■ 



which are not aheady hnked. This is repeated until no unconnected ends remain, in the 
unhkely event that there remain unconnected ends at nodes which are aheady joined, the 
degree of these nodes is reduced. The graph is then 'shuffled' by picking two hnks, g — h, 
i—j at random and changing it to either g — i, h—j ot g—j, h — i with equal probability. 
This is repeated until on average each link has been rewired ~ 10 times. This shuffling 
is intended to remove any bias due to not allowing multiple linkings [2Z|. The P{k)s are 
chosen so that there exists a large connected component (giant component) containing 
almost all nodes, nodes not in the giant component are discarded when calculating the 
magnetisation since the solution methods described above are for the giant component 
only. Network creation took insignificant time in comparison to simulation of the Potts 
model. The Potts model is simulated using the Swendsen-Wang algorithm [2H] ■ System 
sizes of n = 20, 000 50, 000 were easily simulated and were sufficient to allow comparison 
with the solutions described above. 

In figure Q we show the results for magnetisation of the bimodal graph as a 
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Figure 2. Plot of M{T) for the q = 10 Potts model on a random graph with power- 
law degree distribution (see equation (|20|l 'l with J /kg = 1, 7 = 3.5 (rightmost curves) 
and 7 = 4 (leftmost curves), kmin — 2, kmax = 100. Crosses are simulations for 
2 realisations of a network with n — 50, 000 nodes, starting from an initially fully 
magnetised state. Lines are the results of the 2 methods described, the solid curves 
are the method introduced here, the dotted lines are the method of Dorogovtsev et. 
al. The vertical lines are the transition from the metastable unmagnetised state to the 
magnetised state, both methods giving the same value of TmetastaUe- 

function of temperature. This case was chosen specifically because the approximation 
of Dorogovtsev et. al. would perform poorly due to the small degree (2) of half of the 
nodes and the large separation of the degrees (2 : 10). As expected, the method of 
Dorogovtsev et. al. gives an incorrect result for the magnetisation although it's value 
for the transition temperature Tc is surprisingly close to the correct value. 

Figure 121 shows the results for the power law case with exponents 7 = 3.5,4 with 
kmin = 2 and kmax = 100. In practice the expected number of nodes with /c > 60 is 
less than 1 for n = 50, 000 and increasing kmax = 100 does not appreciably change the 
results of our method. Again, the method of Dorogovtsev et. al. differs from the one 
introduced here and from the simulations. In particular, the transition temperature is 
incorrect by ~ 3%. The simulations do not agree quite as well with our solution as 
for the bimodal case, particularly near to the transition temperature. This is expected 
because, as mentioned earlier, the convergence of the iterative procedure slows down 



Potts Model On Random Trees 



8 



close to the transition. A signature of this is that even with such large system sizes, the 
curves obtained for two different realizations of the graph differ slightly. In addition, 
the agreement improves as n increases up to 50, 000. 

Also shown in figures ()1|2|) is the transition from the metastable unmagnetised state 
at a temperature TmetastaUe less than T^. Tmetastabie can be found exactly [21] 

T J/ In (fe^) + - 2) (fc) .21) 

-'-metastable <-'/ i^-"- (fc^) 2 (A;) V / 

For the continuous transitions, i.e. for q = 2, there is no metastable region and so 
Tmetastabie = Tc and thus the method of Dorogovtsev. et. al. gives Tc exactly. 

For these results we have used population sizes of 10^ and, near Tc, 10®. Near 
Tc and Tmetastabie ~ 500 iterations per member of the population were used to ensure 
convergence, less elsewhere. Between Tc and Tmetastabie the magnetisation found depends 
on the initial values of hi, we started from both hi ^ and hi ^ 1, with a slight spread 
in the distribution of hi values. The timescales to produce a full M(T) curve using a 
standard PC were of the order of minutes for the method of Dorogovtsev et. al. and 
hours for the method described here. 

In figures (jH} we show the distributions of the effective fields Xi on neighbouring 
nodes at various temperatures below the phase transition. It is clear that, except in the 
fully magnetised or unmagnetised states, there is a wide spread of local field values. In 
particular, for the bimodal case, the plateau at ?/ ~ 0.83 is because the k = 2 nodes are 
unmagnetised whilst the k = 10 nodes remain magnetised. This is also consistent with 
the form of the M{T) curve in figure (^. 

Figure Q shows the temperature dependence of the relative spread of local 
fields, defined as the ratio between the standard deviation a of the effective fields on 
neighbouring sites hi divided by (h) for (h) > 0. This is non-zero below the phase 
transition, and it attains a maximal value at the phase transition. This shows that 
the approximation of Ref. is far from exact. In spite of this, we found that 

the discrepancy of the results of Ref. [21] with the method described here is almost 
insignificant when the minimal degree kmm is large (say kmin ~ 10). The discrepancy 
between Dorogovtsev et al.'s approximation and numerical simulations also improves 
slowly when q decreases, their result being exact for q = 1 [21] ■ The approximation 
allows some analytic results, e.g. for Tc, to be derived within it and it is numerically 
less demanding and thus significantly faster than ours. 

In conclusion, we have presented a method for solving the Potts model on a random 
graph of given degree distribution. The method does not make the approximation of 
Ref. [25 and thus is accurate in all cases when the local tree approximation is applicable, 
giving results in good agreement with our simulations. 
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Figure 3. The standard deviation of the effective fields on a neighbouring node, 

hi, divided by (h) for (h) > as a function of temperature. Note the drops at the 
transition temperatures, from left to right, the transitions are for: powerlaw 7 = 4 
q = 10 kmin = 2 kmax = 100, bimodal q = 10 k = 2, 10, bimodal q = 10 k = 10, 20, 
and bimodal q = 2 k = 2,10. The insets show the cumulative distribution of effective 
fields X on a neighbouring node, p{x')dx' for the bimodal random graph with q = 2 
and A: = 2, 10 (left inset) and for the powerlaw 7 = 4 g = 10 kmin = 2 kmax = 100 
case (right inset). Various temperatures are shown, T increases from left to right. The 
plateau at y w 0.83 for the bimodal cases is where the k = 2 nodes are unmagnetised, 

Vplateau ~ 

10/(2 + 10). 



References 

[1] M. Newman, in Handbook of Graphs and Networks, edited by S. Bornholdt and H. G. Schuster 

( Wiley- VCH, Berlin, 2003a). 
[2] M. Newman, S. Strogatz, and D. Watts, Phys. Rev. E 64, 026118 (2001). 
[3] F. Ball, D. MoUison, and G. Scalia-Tomba, Annals of Applied Probability 7, 46 (1997). 
[4] D. Watts and S. Strogatz, Nature 393, 440 (1998). 

[5] D. Watts, Small Worlds: The Dynamics of Networks between Order and Randomness (Princeton 

University Press, 2003). 
[6] D. S. Price, J. Amer. Soc. Inform. Sci. 27, 292 (1976). 
[7] A. Barabasi and R. Albert, Science 286, 509 (1999). 
[8] S. Dorogovtsev and J. Mendes, Adv. Phys. 51, 1079 (2002). 
[9] M. Newman, SIAM Review 45, 167 (2003b). 
[10] R. Albert and A. Barabasi, Rev. Mod. Phys. 74, 47 (2002). 



Potts Model On Random Trees 



[11] M. Newman and D. Watts, Phys. Rev. E 60, 7332 (1999). 

[12] D. Callaway, M. Newman, S. Strogatz, and D. Watts, Phys. Rev. Lett. 85, 5468 (2000). 

[13] R. Cohen, K. Erez, D. ben Avraham, and S. Havlin, Phys. Rev. Lett. 85, 4626 (2000). 

[14] R. Pastor-Satorras and A. Vespignani, Phys. Rev. Lett. 86, 3200 (2001). 

[15] R. Kulkarni, E. Almaas, and D. Stroud, cond-mat/9905066 

[16] L. de ArcangeHs and H. Herrmann, Physica A 308, 545 (2002). 

[17] G. Bianconi and M. Marsih, Phys. Rev. E 70, 035105(R) (2004). 

[18] G. Abramson and M. Kuperman, Phys. Rev. E 63, 030901 (2001). 

[19] H. Ebel and S. Bornholdt, Phys. Rev. E 66, 056118 (2002). 

[20] A. Goltsev, S. Dorogovtsev, and J. Mendes, Phys. Rev. E 66, 016104 (2002). 

[21] S. Dorogovtsev, A. Goltsev, and J. Mendes, Eur. Phys. J. B. 38, 177 (2004). 

[22] M. Leone, A. Vazquez, A. Vespignani, and R. Zecchina, Eur. Phys. J. B 28, 191 (2002). 

[23] R. Baxter, J. Phys. C 6, L445 (1973). 

[24] R. Baxter, Exactly Solved Models in Statistical Mechanics (Academic Press, 1982). 
[25] G. Ehrhardt, M. Marsili, and F. Vega-Redondo, inprepar ation] 
[26] R. Cohen and S. Havhn, Phys. Rev. Lett. 90, 058701 (2003). 

[27] R. Milos, N. Kashtan, S. Itzkovitz, M. Newman, and U. Alon, [cond-mat/0312028| 
[28] R. Swendsen and J. Wang, Phys. Rev. Lett. 58, 86 (1987). 



